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Abstract 

We present a new fast solver to calculate fixed-boundary plasma equilibria in toroidally axisymmetric ge- 
ometries. By combining conformal mapping with Fourier and integral equation methods on the unit disk, we 
show that high-order accuracy can be achieved for the solution of the equilibrium equation and its first and 
second derivatives. Smooth arbitrary plasma cross-sections as well as arbitrary pressure and poloidal current 
profiles are used as initial data for the solver. Equilibria with large Shafranov shifts can be computed without 
difficulty. Spectral convergence is demonstrated by comparing the numerical solution with a known exact 
analytic solution. A fusion-relevant example of an equilibrium with a pressure pedestal is also presented. 

Keywords: Grad-Shafranov, plasma physics, Poission solver, spectrally-accurate, conformal mapping, 
high-order, Kerzman-Stein 



1. Introduction 

High performance numerical equilibrium solvers are crucial to the computational effort in plasma physics 
for magnetic fusion applications. Numerically computed plasma equilibria are needed as an input to the 
algorithms responsible for calculating the macroscopic stability and transport properties of plasma config- 
urations. Equilibrium solvers need to be fast, so that computational time can primarily be spent on the 
stability and transport calculations. High accuracy is also required in order to minimize error propagation. 

The equilibrium magnetic configuration in all toroidally axisymmetric magnetic confinement devices, in- 
cluding the tokamak and the spherical tokamak (ST), is determined by solving the Grad-Shafranov (GS) 
equation [TH [251 ISZ] (also known as the Grad-Schliiter-Shafranov equation) . The Grad-Shafanov equation is 
a nonlinear elliptic partial differential equation, which, in general, can only be solved numerically. Because 
of its importance in computational plasma physics, a wide range of numerical methods have been developed 
PP] . These schemes generally fall into two categories. Eulerian or "direct" solvers use a prescribed mesh to 
calculate the unknown function [T31 HZl [121 US US] i while Lagrangian or "inverse" solvers find the mapping 
of the plasma geometry in terms of magnetic coordinates [71 [161 123 EH 112] • The advantages and disadvan- 
tages of one formulation as compared to the other depend on the application of interest [43] , on the plasma 
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geometry, and on the type of inputs that plasma stabihty and transport codes require. For the Eulcrian 
formulations, finite difference [131 HZ] and finite element [T5J [121 [211 I2S1 130] methods are systematically 
favored. 

In this article we present a new, high-order, fixed boundary, direct GS solver and demonstrate its effectiveness 
as part of an Eulerian solver. High-order accuracy is achieved through the combination of three key elements: 
a rescaling of the unknown which reformulates the GS equation as a nonlinear Poisson problem, a spectrally- 
accurate numerical method to compute the conformal map from the plasma cross-section to the unit disk, 
and finally, a fast, high-order Poisson solver on the unit disk. 

Conformal mapping techniques in the context of numerical equilibrium solvers in plasma physics were first 
considered by Goedbloed (TUllIIllH] as a convenient way to decouple the numerical issues associated with the 
plasma geometry from the rest of the problem. We show that for smooth-boundary plasma cross-sections, 
spectrally-accurate conformal maps to the unit disk can be efficiently computed. It is important to note that 
the map needs to be computed only once for a particular geometry, and can be used to iteratively solve the 
nonlinear equation or to compute equilibria for several plasma profiles. By combining this mapping with a 
Green's function-based Poisson solver on the unit disk, we obtain high-order accuracy for the solution as 
well as its derivatives. This is one of the main motivations for this work, since the stability and transport 
properties of the plasma are very sensitive to quantities such as the local magnetic shear and the magnetic 
field curvature which depend on the second derivatives of the GS solution [S] . 

The structure of this article is as follows. In Section 2, we introduce the Grad-Shafranov equation and give 
the mathematical formulation of the fixed-boundary problem. In Section 3, we present a rescaling of the 
unknown which reformulates the GS equation as a nonlinear Poisson problem, and describe the iterative 
method used for its solution. Section 4 contains the details of the numerical conformal mapping scheme 
which maps the plasma domain to the unit disk. The conformal map is computed by solving the Kerzman- 
Stein integral equation, and using its relation to the Szego kernel. The fast Poisson solver on the unit disk 
is presented in Section 5. In Section 6, we illustrate the efficiency and accuracy of our new solver with a 
few specific examples, and in Section 7 we conclude with a discussion of the limitations of the solver, and its 
potential extensions. 

2. The Grad-Shafranov equation as a nonlinear eigenvalue problem 
2.1. The Grad-Shafranov equation 

In toroidally axisymmetric systems, in the usual cylindrical coordinate system (R,(j),Z), the magnetic field 
can be expressed as 



where is a unit vector in the toroidal direction. By the assumption of axisymmetry, none of the physical 
functions of interest depend on the angle ip. Here 2TT'i>{R,Z) represents the poloidal magnetic flux, and 
27rg(5') = — 7^(5') is the net poloidal current flowing in the plasma and the toroidal held coils. The flux 
function 'J satisfies the Grad-Shafranov equation 



B = 




(1) 




(2) 



2 



where p is the plasma pressure. Both p and g are appHcation-specific functions of 5*, which, along with the 
boundary conditions, determine the equilibrium. Indeed, once ([2| is solved and 5* is known, p and g can be 
immediately evaluated, B can be computed from , and the current density J in the plasma is given by 
(see [5, for example): 

2.2. Boundary conditions 

Equation ([2| is a second-order elliptic nonlinear partial differential equation. Depending on the boundary 
conditions to be enforced, we distinguish between two general classes of problems: (i) fixed-boundary prob- 
lems in which the plasma boundary is prescribed, with 5* = constant on the boundary and one solves for 
\l/ inside the plasma, and (ii) free-boundary problems where the current flowing in a set of external coils 
is given and one has to find ^' such that the equilibrium is self-consistent with these currents. In this pa- 
per, we consider only fixed-boundary problems. That is, the shape of the plasma boundary is given and 
^! = constant on the boundary. Furthermore, since only derivatives of ^ represent physical quantities, we 
can choose ^ = on the boundary without loss of generality. 

Our motivation for focusing on the fixed-boundary problem is two-fold. First, a large number of plasma 
stability and transport numerical codes take fixed boundary equilibria as their initial data. In particular, 
parametric studies to understand and optimize plasma properties by varying a few parameters defining a 
given generic plasma boundary are very common. Second, many free-boundary GS solvers use iterative 
schemes which require a robust fixed-boundary solver in the iterative loop |20l I40j . 



2.3. The eigenvalue problem 

As explained in the previous section, the generic form of the fixed-boundary GS equation can be written as 

(4) 

* = on dfl, 

with 

where n represents the interior of the plasma, and dfl the plasma boundary. The nature of the previous 
problem depends on the behavior of the right hand side as a function of ^I^; the user-specified plasma profiles 
determine this. For many profiles, the right hand side F is of the form 

= (5) 

In this case, the function ^I* = is a trivial, but not physically relevant, solution to equation Q. It is 
well-known [TTJ [511 SO] that in this case the Grad-Shafranov equation is solved as an eigenvalue problem. 
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which is hnear when F is linear in and nonlinear otherwise. In particular, consider the scalings 

' ' I^U d* ' 

d5^ ^d5^ Fr* i? = — F(* i? Z) (7) 

where \^\m = sup|^'| and cr is a normalization factor for the pressure and poloidal current profiles. Note 
that 'I', as defined in ([6|, takes values in the interval [—1, 1]. Inserting the normalizations in equations ([6| 
and into Q, we find 

Equation ([s]) is in the desired form, highlighting a well-known property of the GS equation: its scale- 
invariance under the transformation 

cr, i?, Z) ^ (A*, A^CT, i?, Z) . (9) 

This scale invariance implies that by defining the ratio a = (t/|'I'|^, we can write the fixed-boundary GS 
equation in the following scale-independent form: 

A*v& = aF(^,R,Z) in ft 

(10) 

^ = on dfl. 

This is now an eigenvalue problem to be solved for the eigenfunction 5* with eigenvalue a. When transforming 
back to the original variables that represent physical quantities of interest, there are two choices. Either 
I^'Ito is given (or alternatively I^, the total toroidal current flowing in the plasma), and the magnitude of 
the pressure and current profiles is computed from the relation a = |*|j^(t, or the normalization a is given 
and the corresponding \'^\m or is determined. 

On the other hand, the pressure profile or the poloidal current profile may be specified to include terms 
which are linear in *. Under this assumption, their derivatives can be written in the form 

-Mo^ = -Mo5(vI/)vl/ + C, -^^Ti^)^ + A, (11) 

with A and C constants such that C + A 0. This choice of p and g result in the right hand side F of 
equation Q having a constant term. This choice corresponds to a discontinous toroidal current at the 
plasma edge, and equation Q does not correspond to an eigenvalue problem. While this is not a physically 
relevant choice, it plays an important role in the benchmarking of numerical Grad-Shafranov solvers. When 
the functions S and T are identically zero (corresponding to the so-called Solov'ev profiles |39| ) exact analytic 
solutions to the GS equation with fixed boundary conditions can be constructed O El |44] , which can then 
be used to test the accuracy of a numerical solver. In Section |6] we present convergence results of our solver 
benchmarked against such analytic solutions. 

In what follows, we describe a new numerical scheme which solves equation Q with high accuracy. Both the 
eigenvalue and non-eigenvalue problems will be considered. For the eigenvalue problem, we only treat the 
question of finding the smallest eigenvalue a and its associated eigenfunction because almost all magnetic 
fusion confinement problems of interest have a single extremum of * within the plasma region. Our procedure 
could however be easily generalized to find other eigenvalues, for so-called doublet or more ornate plasma 
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configurations. 



3. The Grad-Shafranov equation as an iterative Poisson problem 

3.1. Transformation to a two-dimensional nonlinear Poisson problem 

We consider the generic form of the fixed-boundary Grad-Shafranov equation, given by equation Q. In aU 
but the simplest cases, an iterative scheme is necessary since either , R, Z) is a nonlinear function of 
or a nonlinear eigenvalue problem is to be solved (or both). As defined earlier, the differential operator A* 
is a second-order elliptic operator, similar to the Laplacian: 

If iteration is to be used, one might consider rewriting equation Q as 

'fBR + '^zz = F{'f,R,Z) + ^'fR inn 

R (12) 

f = on dCl, 

and solve iteratively using a two-dimensional Cartesian Poisson solver at each step, holding the right hand 
side fixed at the previous iterate. However, this approach has one caveat. It is empirically observed that 
fixed-point iteration converges faster when the right-hand side is slowly varying. It would therefore be 
preferable to have a right-hand side that is only dependent on the unknown function itself, and not on its 
derivatives. In order to obtain this desired form, we consider the following scaling of the unknown function 

UiR,Z) = ^^iR,Z). (13) 
This transformation is clearly singular at i? = 0. Fortunately, this singularity does not have any consequence 



because in all applications of interest, the physical domain excludes the Z-axis. Inserting transformation ( 13 ) 
into Q, we readily find that U satisfies the following equation: 

Urr + Uzz =HU,R,Z) inO 

(14) 

[/ = on ar2 , 



where J- is defined as 



T{U, R, Z) = -^F{VRU, R, Z) + -^U. 



Here, is a function of U only, and not of its derivatives. Equation ( 14) has the form of a two-dimensional 
nonlinear Poisson problem (treating R and Z as Cartesian coordinates), which we solve iteratively as de- 
scribed in the next section. Once U and all its derivatives have been computed, all quantities of physical 
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interest can be evaluated to the same accuracy using the relations: 

ZV -ft 

"^zz = ^Uzz- 



3.2. Iterative solution for the Grad-Shafranov equation 



When the nonlinear Poisson problem defined by (14) does not correspond to an eigenvalue problem (i.e 



the right hand side contains a constant term), we use a straight-forward fixed-point iteration scheme. At 
each iteration, the right-hand side of ( [l4| is evaluated at the previous value of U, and the Laplacian A = 
d^/dR^ + d"^ jdZ'^ is inverted to calculate the updated value of U . In other words, on the V'^ iteration, we 
solve 

AC/j =^([/,_i,i?,Z) (15) 

for Ui. The iterative process is stopped when \\Ui — C/i-i||oo < e for a few consecutive iterations, with e small 
and chosen by the user. 



The iterative scheme for the related eigenvalue problem is more subtle. Equation ( 10 ) describes a nonlinear 
eigenvalue problem to be solved for the smallest eigenvalue a and its associated eigenfunction ^, with 
ll^lloo = 1- Here, the norm || • ||oo has its usual meaning, i.e. ||^'||oo — supd^fj). When solving a linear 
eigenvalue problem, the eigenvalue does not depend on the scaling of the eigenfunction. However, this is 
generally not true for nonlinear problems. Therefore, we use a modified version of the well-known inverse 
iteration method |41j that evolves the eigenvalue-eigenfunction pair simultaneously instead of evolving the 
eigenfunction alone. First, the previous iterate eigen-pair (^f;-!, (T;_i) is used to solve the PDE 

A**z =a,_iF(§z_i,i?,Z) inn 

(16) 

= on dfl . 

Then, the next eigen-pair iterate {"^i^ai) is computed via correcting by the norm of 

Once again, the iterative process can be terminated when ||\['/ — 5';_i||oo < e for a few consecutive iterations, 
with e small and chosen by the user. 



Applying the transformation (131 to (16) and (17) gives the iterative scheme in terms of U. First we solve 
the Poisson problem for Ui : 

AUi = ^FiVRUi-i,R,Z) + -^Ui-i in ft 
Ui = on dn . 
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Figure 1; Conformal map of a domain to the unit disk. 



Then, we compute the updated Ui and ai as: 



U, 



RUi\ 



o-i-i 



RUi\ 



(19) 



This eigenvalue iteration scheme generally converges at a geometric rate, as shown in [35 . Note that in this 
scheme, we calculate lI'v/RWfUoo by first finding the maximum on the computational grid, and then using 
Newton's method to refine the maximum. Without this latter step, the problem we solve would depend on 
the computational grid and would converge only linearly, and not spectrally approximate ||\/i?C/i||oo as the 
grid is refined. 

In summary, we have shown that the Grad-Shafranov equation can be solved using a sequence of iterations 
which require a single two-dimensional Poisson boundary value solve in Cartesian coordinates at each step. 
We now move on to developing a highly accurate Poisson solver in with Dirichlet conditions [/ = on dQ,. 
This will be accomplished in two steps. First, we present a fast spectrally-accurate method for numerically 
computing the conformal map from our domain Q, to the unit disk. Then, we develop a fast, high-order 
Poisson solver on the unit disk using separation of variables. 



4. Spectrally-accurate conformal mapping for smooth boundaries 

Computation of the forward map of the boundary to the unit circle 

In this section, the notation will be as follows. The original (i?, Z) domain will be identified with the complex 
z-plane, with a point in this plane denoted as z — x + iy. The image domain will be referred to as the w- 
plane, with a point denoted as w ^ a + ij3. The forward map is the complex valued function w = W{z) 
with real and imaginary parts a — A[x^ y) and /3 = B{x, y). The inverse map is the complex valued function 
z — Z{w) with real and imaginary parts x — X{a, /3) and y = Y{a, (3) (see Figurejl]). Note that the function 
Z is obviously not the same mathematical object as the coordinate Z of the (i?, ((>, Z) coordinate system 
used in this article; this notation should not lead to any confusion since the two objects are never used 
simultaneously and can easily be distinguished based on the context. 

The Riemann mapping theorem guarantees the existence and uniqueness of an analytic map between any 
connected set and the unit disk. It is straightforward to show that under such a mapping, solving the Poisson 



7 



problem 



Au{x,y) = f{x,y), u\ 



= 



in the original domain $7 is equivalent to solving the scaled Poisson problem 



on the mapped domain Di (in our case the unit disk), with 



dZ 



dw 



= 



u{x, y) = v{A{x, y),B{x, y)). 

We compute the forward map W on dfl using the Kerzman-Stein integral equation, see [22l[23j. In particular, 
the derivative of the forward conformal map W on dfl is given by 



W'{z) = 2tt 



S{a, a) 



S{a, a) 



(20) 



S{z, a)S{z, a)ds2 



an 



where a is the point mapped to the origin, i.e. W{a) = 0. In the examples that follow, we always choose a 
as the geometric center of fl (see Figure [6]). The function S is the Szego kernel, which satisfies the integral 
equation 

S{z,a)+ Aiz,t)S{t,a) dst = H{a,z). (21) 



J so 

Here, the overline represents the complex conjugate, H is the Cauchy kernel, 

, \ dz/dl 
H{w,z) = — , 

ZiTl Z — W 

and the kernel A is known as the Kerzman-Stein kernel, defined as 



A[z,w) 



H{w, z) — H{z,w), z^w 
0, z — w . 



Numerical implementation of (21) is straightforward since the kernel A is smooth. Since the trapezoidal 
rule is spectrally-accurate when applied to sufficiently smooth periodic functions, the boundary 917, and 
therefore the integral in (21), can be discretized using A^i points equally spaced in arc- length. The number 
of discretization points A^i is chosen so that the boundary parameterization and its tangential derivatives 
are fully resolved to machine precision. Letting Zk and {dz/dl)k denote equispaced points and tangential 
derivatives on dfl (specified either directly, or obtained from an equation of the plasma surface) , the integral 
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equation in (21 ) is then discrctizcd resulting in a linear system of the form Ms = b, with 




1, 



h l{dz/dl)j 
27ri Zj — Zk 



+ 



{dz/dl)k 



J 



1 \{dz/dl)j 



27ri Zj — a 
S{zk,a), 



where h = L/Ni is the step size (with L denoting the length of the boundary d^l). Naive inversion of the 
linear system Ms = b to obtain s requires 0(-/Vf ) operations. There exist fast methods ^34 which can be 
used to accelerate the solution of this linear system. However, since the conformal map needs to be computed 
only once in our GS solver, and since the iterative PDE solver is a significantly more expensive step, we have 
relied in our experiments on the simpler dense matrix formulation. 

Once s is known, the derivative of the forward map for the points on the boundary can be calculated 



Since s is known at equispaced points in arc length, the FFT can be used to compute W{zk) by Fourier 
integration. The integration constant is determined so that the unit disk is centered at the origin. Specifically, 
three points on the mapped circle can be used to find its center, which is then subtracted from the map. 
Note that by using Fourier differentiation, we can also compute W" , which is necessary for the mapping of 
the second derivatives of functions defined on the original domain and the mapped domain |35j . 

An example of the boundary values of the forward conformal map is shown in Figure |2] It is easy to see that 
for arbitrary boundaries, equally spaced points on the boundary dfl do not map to equally spaced points on 
the unit circle, a well known feature of conformal mapping known as crowding. In fact, the level of crowding 
is exponential in the aspect ratio of |34j . In many cases, high-order Poisson solvers on the disk which take 
advantage of Fourier methods require that discretization points be equally spaced in the angular direction. 
Since these points do not coincide with the ones obtained via the forward conformal map, it is necessary to 
resample the mapped points. 

4-2. Computing the inverse map for equally spaced points on the unit circle 

To obtain a uniformly spaced angular grid on the unit disk, we proceed in two steps. First, the boundary 
values of the forward map : — > Di are oversampled via the FFT. This is shown in Figure [3j The 
amount of oversampling necessary is at least the crowding factor of the map. 

Second, the image points on dDi (the unit circle) are resampled to equispaced ones, effectively inverting 
the conformal map. We consider z e dfl as a function of w G dDi and use Lagrange interpolation to find z 
values corresponding to equally spaced w values on the circle. In particular, barycentric interpolation is used, 
which behaves well even for high order [3]. For example, 8*''-order barycentric interpolation uses as data 
4 points on the left and the right of the target interpolation point. The oversampling rate in the previous 



using (20): 





(22) 
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Figure 4: Inverse of the conformal map of the boundary — equally spaced on the circle: original domain (left), unit circle 
(right). 

step can be tuned so that this interpolation procedure gives the desired accuracy, sometimes requiring a rate 
larger than the crowding factor. The end result of the interpolation step is shown in Figure |4] 

The equally spaced w boundary points we obtain in this step will be used for the angular grid of a high- 
order Poisson solver. We know their image z = Z(w) on the boundary, but do not yet know the image 
of interior points of the grid. To compute the inverse map Z (as well as Z' and Z") for interior points, 
the Cauchy integral formula can be used. Since the boundary of the domain is the unit circle, the integral 
formula is particularly easy to evaluate with the FFT. The inverse map of the interior points of a regularly 
spaced Nj. x Ng grid on the unit disk can be computed in ©(iV^A^e log A^e) time, the same computational 
complexity as the Poisson solver on the disk (as shown in the next section). This method is spectrally- 
accurate, and works all the way up to the boundary without any loss of precision or any requirement for 
adaptive integration (due to the proximity of the singularity of the Cauchy integral to the contour, see [3S] 
for a detailed discussion). 

We end this section with a short discussion on the previously mentioned issue of crowding, inherent to the 
conformal mapping technique, and clearly visible in Figures [2] [sj and|4] (see also [SS])- For certain shapes 
(such as ones with a high aspect ratio), conformal maps are clearly not a viable method, at least when 
used in combination with solvers which require an equally spaced grid on the disk. For example, a 10 to 
1 aspect ratio ellipse has a crowding factor on the order of 10^", which would obviously lead to prohibitive 
oversampling. Fortunately, for the cross-sections of tokamaks and spherical tokamaks, the crowding factor 
is much more manageable, on the order of ^ 10 - 20. For numerous geometries relevant to magnetic fusion, 
conformal mapping is an efficient method, as will be shown in more detail in Section [6] 

5. High-order Poisson solver on the unit disk 

What remains to be discussed in order to complete our Grad-Shafranov solver is the solution of the two- 
dimensional Poisson equation on the unit disk. There are, of course, a number of fast solvers available for 
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this problem. For the sake of completeness, we describe one here that is direct, high-order accurate, and 
straightforward to implement. Other high-order schemes can be found, for example, in [51 138] . 



5.1. Separation of variables and boundary conditions 
Let us consider the fixed-boundary Dirichlet Poisson problem 

Aw = / in Di 

(23) 

V = Q on dDi , 

where Z)i is the unit disk. Using separation of variables in the usual polar coordinates {r,0), we represent 
both the solution u and the right-hand / side as Fourier series: 

oo oo 

v{r,e)= ^ ^)„(r)e»^ f{r,9)= ^ /„(r)e™^ (24) 



Substituting these expressions into ( 23 1 , we have the following radial ordinary differential equation for each 
Fourier mode n: 



1 



+ ^v'nir) - '\vn{r) = fn{r) 



r 



(25) 



*„(!) = 0. 



In order for equation (25) to be well-posed, a second boundary condition must be enforced. We obtain this 



condition by requiring regularity of the solution as r — )■ 0. For the n = mode, multiplying (25) by r and 
taking the limit r — >■ leads to the condition 'C„(0) = 0, under the assumption that /„, f)„, and {)„ are 
bounded. For n 7^ 0, multiplying the equation by and again taking the limit r — leads to ^^(O) = 0, 
under the same boundedness assumptions. In summary, we are to solve the following decoupled system of 
equations in the radial direction: 

^nir) + -v^{r) :^Vn{r) = fn{r) , z;„(0)=0, {)„(1) = for n = , 

\ 1 (26) 

// 1 / Tl 

+ ~ ^Vn{r) = fn{r) . ^)„(0) - , i)„(l) = torn ^ . 
5.2. Green's functions solution to the radial equation 



Equation (26) is a well studied ODE in mathematical physics, whose Green's function is known. Using 
convolution with the Green's function, a particular solution to the ODE which does not satisfy the 
boundary condition at r = 1 can immediately be written down. Then, a correction -C^ can be found which 
solves the homogeneous equation, and such that the sum Vn = + Vn satisfies the ODE and the boundary 
conditions, both at r = and at r = 1. To this end, we proceed as follows. 
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The Green's function with the proper behavior at and infinity for the ODE in (26) is 



{s log s , r < s 
n = 0, 
slogr , r > s 



2\n 

and the convolution solutions are 



■r"l"lsl"l+i , r>s 



CXD 



t)^(r)=logr/ sfn{s)ds+ s\ogsfn{s)ds n = 0, 

Jo Jr 

Since we are solving the Poisson problem on the unit disk, we can set /„(?") = for r > 1, and the convolution 
solutions to the inhomogeneous equation can be written in the form 



v^{r) = \ogrj sfn{s)ds + j s\ogsfn{s)ds n = 0, (27a) 

(^r-l"l^''sl"l+7„(s)ds + rl"l^'s-|"l+V„(s)ds^ n^O. (27b) 



The solution to the homogeneous equation 

1 . 

inir) + -Vnir) - ^vnir) ^ (28) 
satisfying the regularity condition at r = is 

^,f(r) = c„rl"l, 

where c„ is a constant to be determined from the boundary condition at r = 1: 



Therefore, the general solution to (26) satisfying the boundary conditions at r = and r = 1 is 

z)„(r) = ^^(r)-t),^(l)rH. (29) 

As previously discussed, one of the major advantages of Green's function methods is that the radial and 
angular derivatives of the solution v can be calculated explicitly by differentiating the separation of variables 
representation. Numerical differentiation is never required. This allows for accurate computation of the first 
and second derivatives of ^E", which are required for stability and transport calculations. 



The radial derivatives are computed by direct differentiation of (27a 27b) and (29). For the n = mode. 
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we find 



i'nir) = fn{r) ^ / s/„(s) ds 



(30) 



For non-zero n modes, we have 
1 



v'ir) 



-\n\-l 



\n\-l 







-|n| + l 



/„(s)ds-H«^(l)rH-i, 



i:ir) = fn{r)- 



\^\ + l„-|n|-2 



2 



fn{s)ds 



(31) 



Angular derivatives are obtained via differentiating tlie Fourier series representation (24 1, i.e. by a simple 
multiplication of in. The first and second partial derivatives of the solution v are then computed according 
to 

OO OO 

71— — OO n— — oo 

OO OO 

n— — oo n— — OO 

OO 



n— — OO 



where the radial derivatives v'^{r) and {""(r) are given in equations (30) and (31). As is well known, we 
observe that taking one derivative introduces a condition number of 0(n), and taking a second derivative 
introduces a condition number of ©(n^), due to the multiplication by n and ■n? respectively. 



5.3. Numerical considerations 
5.3.1. Grid setup 

Since the Poisson solver previously outlined relies on separation of variables in polar coordinates, we build 
a tensor grid in the radial and angular variables r and 6. In the 9 variable, we compute and sum Fourier 
series using the FFT. When sampled at equispaced points, the representation in terms of Fourier series is 
spectrally-accurate for smooth data. Let Nq denote the number of equispaced angular grid points in 9. In the 
radial direction, however, the domain is not periodic. Therefore, an appropriate high-order representation is 
a piecewise Chebyshev grid: the interval [0,1] is divided into Nl intervals, and a P*'' order Chebyshev grid is 
constructed on each interval. The total number of points in the radial direction is therefore Nj. — PN]^. The 
order P can be chosen as desired; in all subsequent examples we set P = 16, yielding a 16"^ order scheme. 
The grid used for a typical tokamak geometry (discussed in more detail in Section |6| is shown in Figure [5] 
along with its image under the inverse conformal map in the original domain. 
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A3z 



A3w 




5.3.2. Evaluation of Green's function convolutions 



Naive implementations of formulas (27a 27b) require 0{N^) work. However, the first integrals in (27a 27b) 
can be computed recursively from r = to r = 1, while the second integrals in (27a 27b) can be computed 
recursively from r = 1 to r = 0. Thus, all integrals can be calculated in 0{Nr) work. 



Furthermore, some care must be taken in the computation of the integrals in ( 27a 27b ) near r = because of 
the rapid growth/decay of the monomials s'"' and s"'"' for large mode numbers n. The difhculties associated 
with floating-point over low/underflow are easily avoided by rescaling. A second problem is that for n large, 
the integrands rl'^l/s*^!"'"^^ and sl"l+^/rl"l /n(s) are poorly resolved by the composite Chebyshev grid. 

Since /n(s) is well-resolved, however, the change of variables cr = (r/s)'"' and fi = (r/s)-'^~l"l for the first and 
second cases, respectively, yield well-resolved functions in the transformed variable. These mapped integrals 
can be computed accurately with standard quadrature rules [35) . In our experiments below, we used a 16"^ 
order Gaussian rule. 



5.3.3. Convergence and run time 

If the data given is smooth, the Fourier series representation is spectrally-accurate. The overall order of 
convergence of the algorithm is therefore P, the order of the piecewise Chebyshev polynomials used in the 
radial ODE solver. 

The run time complexity of the algorithm is 0{NrNg logiVg), nearly optimal with respect to the number of 
grid points. In detail, we compute 0{Nj.) FFTs of size Ng at a cost of 0{Nj.Ng \ogNg), and Ng ODE solves 
of complexity 0{Nr) are performed at a total cost 0{NrNg). Note that here we treat P as a fixed constant 
that would not be increased if the grid is refined. 
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R 
Ro 



R /R^ 

max 



Figure 6: Geometric definition of the parameters e, K, and 5. 



6. Numerical tests Examples 



Ro 



6.1. Comparison with exact analytic solutions 

In order to test the accuracy of our Grad-Shafranov solver, we first consider a case where exact solutions are 
known. As discussed in Section [2] these are easy to construct for profiles of the form given in equation ( 11 ) 
with 5 = and T = 0. For simplicity we set A = Q (corresponding to a plasma which is neither paramagnetic 
nor diamagnetic [5]) and normalize the pressure such that C = 1. These choices are equivalent to solving 
the GS equation 

A*«' = i?2^ (32) 



which we solve both numerically and analytically. A simple analytic solution to ( 32 ) which is relevant to 



magnetic fusion can be constructed by following the methodology given in [5] . We sum a particular solution 
to the equation, i?^/8, with three solutions to the homogeneous equation: 



8 



(33) 



The free coefficients di, d2, and d^ are determined so that the contour ^E* = represents a reasonable 
plasma cross-section. Specifically, we introduce three characteristic quantities describing the cross-section of 
a magnetic confinement device, as shown in Figure [6} the inverse aspect ratio e, the elongation k, and the 
triangularity 6. 
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Solov'ev sohition up-down symmetric (e=0.32, k=1.7. (5=0.33) 

0.6 
0.4 
0.2 


-0.2 
-0.4 
-0.6 

-1 -0.5 0.r> 1 

Figure 7: Solution to (j32| - ITER parameters. 

The boundary conditions to be enforced are: 

^'(1 + e, 0) = Condition at the outboard midplane, 
^'(1 — e, 0) = Condition at the inboard midplane, 
^'(1 — Se, ne) — Condition at the top/bottom, 

which gives the following system of three equations for di, 3,2 and ^3: 



I 
I 



0.005 

0.01 

0.015 

0.02 

025 

0.03 

0.035 



(34) 



1 


(n 




(1 + ^)* 




" di " 


1 

^ ^8 


" (1 + e)^ " 


1 


(1- 




(1-6)4 






(1-6)4 


1 


(1- 


5ef (1- 


-5e)4-4(l-5e)2K2e2 _ 




. da . 




. (1 - Se)^ _ 



(35) 



Equation (35) is easily inverted, and once the coefficients di, ^2, and ^3 are determined, the analytic solution 



to (32 1 given by (33) is straightforwardly computed. Furthermore, the boundary of the plasma is given by 



the equation 



i?4 



rfi + rfai?^ + dsiR'^ - AR^Z^) = 0. 



Using this plasma boundary allows the numerically obtained solution to be compared with the exact analytic 
one. The grid can then be refined enabling us to check the convergence of the scheme. Two examples are 
shown here: an ITER-like case [HE] with parameters e = 0.32, k = 1.7, S = 0.33, and an NSTX-like case 
[55] with 6 = 0.78, K = 2, S = 0.35. Convergence is measured by computing the L°° norm (i.e. sup norm) 
of the difference between the numerical solution and the exact solution. For the ITER-like case, contour 
plots of the numerically obtained solution are shown in Figure [7| and the convergence behavioris shown in 
Figure |8j For the NSTX-like case, the analogous figures are [9] and [TO] 

Figures [Sj and 10 demonstrate the spectral convergence of the solution and its first and second derivatives. In 
the ITER-like case, ^E* is computed with an accuracy close to machine precision for a 600 by 600 grid. One can 
notice that the errors are a bit larger for the first and second derivatives. This is because we solve the Grad- 
Shafranov equation with Dirichlet boundary conditions, and in order to satisfy the differential equation near 
the boundary we effectively compute two derivatives of the boundary condition. This procedure introduces 
numerical differentiation errors, albeit only of size 0{n) and 0{n'^). 
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Convergence for the np-down symmetric Solov'cv ease (e— 0.78. 2.0, (5—0.35) 
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Figure 10: Convergence of the numerical solution to ((321, and of its first and second derivatives - NSTX parameters. 



One can also observe that the NSTX-hke case has worse convergence than the ITER-like case, and requires 
a finer grid to obtain comparable precision. This is a direct consequence of the conformal mapping method. 
Indeed, the NSTX-like cross-section has a crowding factor of about 20, compared to about 7 for the ITER- 
like case. Therefore, an extra oversampling factor of 3 is required in constructing the conformal map. In 
addition, the NSTX-like boundary requires slightly more points to be resolved. 

The GS solver described is obviously not limited to problems with up-down symmetry (symmetry about 



the line Z = 0). Similar numerical tests were done with exact up-down asymmetric solutions to (32 1, and 
showed very similar performance |35| . Not suprisingly, the solver is also not limited to low pressure cases 
(low-/?). Equilibria with significant Shafranov shifts have been computed without any difficulty |35j . 



6.2. Numerical equilibrium with pressure pedestal 

As an illustration of the Grad-Shafranov equation as a nonlinear eigenvalue problem (see Section [2]) , we 
consider a generic pressure profile corresponding to an equilibrium with a pressure pedestal [H [ST] : 

p(vl/) = (Ci + (l - e-*'M , 



where Ci and C2 are normalization constants, set to Ci = 0.8 and C2 = 0.2 in the examples shown below. 
The parameter 77 is a constant associated with the width and the steepness of the pressure pedestal. Equilibria 
are calculated for three different pedestal steepnesses corresponding to 77 = 0.1, 77 ~ 0.02 and i] = 0.005 (see 
Figure 111. Assuming, as before, that the plasma is neither paramagnetic nor diamagnetic {dg^/<N> = 0), 



the GS equation becomes the following nonlinear eigenvalue problem (to within a constant factor that is 
subsumed into the eigenvalue): 



2C2 



* (1 - e-*'/") + -(Ci + C72*')*e-*'/'' 



(36) 



We solve ( 36 1 for the eigenfunction \1/ and the eigenvalue a using the iterative algorithm presented in Section|3] 
The contours of the solution for the case 77 = 0.005 are shown in Figure [T2| 
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Choices of preKsiire profiles Clioices of dp/d'^ 




Figure 11: Pressure profiles for rj = 0.005, 0.02, 0.1. 
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In order to understand the influence of the steepness of the pressure pedestal on the solution, 5" and its first 
and second radial derivatives are plotted along the line z = 0. The results are displayed in Figure [TSj We 
see that the solution ^ and its first derivative depend only weakly on rj. However, in the vicinity of the 
edge of the plasma, in the pedestal region, the behavior of the second derivative '^/rr is very sensitive to the 
value of rj. 

7. Conclusion 

This article describes a new, fixed-boundary, direct Grad-Shafranov solver that relies upon conformal map- 
ping and Green's function methods to compute high-order accurate plasma equilibria. Its attractive features 
are its proven spectral convergence, its speed, and its versatility in that arbitrary plasma boundaries can be 
given as input (both up-down symmetric and asymmetric). The solver is shown to have very good perfor- 
mance for a wide range of pressure profiles, including fusion-relevant profiles with a steep pressure pedestal. 
Equilibria with large Shafranov shifts are also computed without difficulty. 

The high accuracy achieved for the first and second derivatives of the solution suggests that this solver could 
succesfully be combined with macroscopic stability codes, or implemented in transport studies on long time 
scales in which the equilibrium profiles are evolved self-consistently according to the GS equation. This 
is particularly true for low aspect ratio devices such as the tokamak, for which crowding in the conformal 
mapping is not very pronounced, and relatively coarse grids yield very accurate solutions. For high aspect 
ratio devices such as the ST, the high crowding factor observed in the conformal mapping suggests that 
alternative methods may be more efficient. In this regime, a Poisson solver relying on fast multipole methods 
[5J [121 mi 133] or compression |31j, both of which are highly adaptive, may represent an attractive alternative. 
This would also permit calculations on non-smooth plasma boundaries with fusion-relevant X-points to be 
performed. These alternative approaches are subjects of ongoing research. 
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